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1.  INTRODUCTION 
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This  report  concludes  the  second  phase  of  a  three-phase  collaborative  effort 
between  the  U.S.  Army  Engineer  Topographic  Laboratories and  the  National 
Bureau  of  Standards* (NBS).  This  phase  is  the  object  of  Interagency  Agreement 
E8786K041  covering  the  time  period  from  January  to  September  1986.  The  first 
phase  of  this  effort  took  place  from  May  to  December  1985  under  Interagency 

ment  E8735K137.  The  third  phase  of  this  effort  is  in  the  planning  stage. 


^The  main  objective  of  this  effort  is  to  develop  and  implement  a  *'contour-to- 
grid  algorithm,  that  is,  an  algorithm  capable  of  converting  digitized  contour 
data  into  Digital  Terrain  Elevation  Data  (DTED).  The  approach  is  to  create  a 
synthetic  terrain  surface  based  on  digitized  contour  information  and  then  to 
calculate  any  grid  from  that  surface.  Five  tasks  are  to  be  performed: 


0  Edit  a  digitized  contour  data  set  (input) j 
Ol)Thin  (reduce,  sample)  the  input  data^ 

<i  ^Triangulate  the  selected  points^ 

Ik) Generate  a  smooth  synthetic  surface^ 
q  J>  5) Generate  the  desired  grid.^ 


-7*- 


£  The  first  phase,  which  was  completed  in  December  1985,  addressed  the  first  three 
tasks. ^  A  description  of -the  work  done  and  the  results  obtained  appear  in 
Witzgall,  Bernal  and  Mandel  (85).  jThe  second  phase,  to  which  this  report 
refers,  elaborates  on  the  previous  results  and  includes  work  in  all  five  areas. 
The  effort  demonstrates  the  feasibility  of__performing  such  tasks  for  large  data 
sets,  such  as  a  Digital  Graphic  Recorder ^data  tape  of  roughly  500,000 
points.  (— - - 


Details  about  the  data  structure  used,  the  problems  encountered,  and  the 
necessary  transformations  will  be  discussed  in  the  second  section  of  the  report. 
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In  the  third  section,  we  will  describe  a  tolerance  band  algorithm  for  thinning 
the  data  set  to  a  size  suitable  for  triangulation  and  subsequent  interpolation. 
This  tolerance  band  algorithm  is  a  one-pass  version  of  the  method  by  Reumann  and 
Wickam  (74)  in  analogy  to  techniques  proposed  by  Williams  (78,81).  It  was  used 
to  create  sets  of  sample  contour  points,  ranging  from  roughly  40,000  to  70,000 
points.  An  algorithm  previously  developed  by  MBS  was  then  employed  to  determine 
the  Voronoi  triangulations  of  these  sets  of  points  in  a  plane  as  illustrated  in 
Figure  1.  This  process  required,  however,  the  development  of  a  decomposition 
algorithm  for  data  sets  consisting  of  more  than  50,000  points.  This  and  other 
necessary  adaptations  are  discussed  in  the  fourth  section  of  this  report. 

Section  Five  will  include  a  detailed  dircusslon  of  the  methodology  employed  to 
generate  the  surface,  including  the  algorithm  developed  to  produce  rectangular 
grids  for  any  given  unit  sizes.  The  testing  procedure,  along  with  a  description 
of  the  structural  layout  used  while  performing  the  different  tasks  appears  in 
the  sixth  section.  This  section  will  also  contain  a  discussion  of  additional 
issues  that  are  expected  to  be  relevant  to  further  work  in  this  area. 

Plots  were  obtained  for  the  contour  lines  in  their  original  form  as  well  as  in 
generalized  form  corresponding  to  40,000  and  70,000  point  samples  using  a  Gerber 
Plotter,  The  Voronoi  triangulation  of  a  40,000  point  sample  was  also  plotted. 
The  surface  generation  software  which  wac  prepared  as  part  of  this  effort,  was 
applied  to  a  a  40,000  point  sample.  The  resulting  surface  was  tested  using  the 
original  data,  independently  of  the  sample  data  set. 

All  computations  were  initially  carried  out  on  a  VAX  11/750  system  and  later 
(July  1986)  transported  to  a  VAX  11/780  system  at  ETL.  They  were  found  to  be  in 
the  range  of  extensive,  but  routine  computations  of  the  kind  that  can  be 
expected  as  part  of  the  normal  load  of  such  systems. 

Triangulation  based  surface  modelling  has  been  considered  for  some  time,  e.g. 
Peucker  and  Douglas  (75).  For  other  rel  fed  work  see  Davis,  Downing,  and 
Zoraster  (82),  and  Grotzinger,  Danielson,  Caldwell,  and  Mandel  (84). 


2.  DIGITAL  DATA  STRUCTURE 


The  first  task,  performed  in  this  effort  consisted  of  obtaining  and  editing  a 
digital  contour  data  set  making  it  compatible  with  the  software  developed  for 
the  other  tasks*  All  the  necessary  computations  were  performed  originally  in  a 
VAX  11/750  system  end  then  transported  to  a  VAX  11/780  with  a  VMS  operating 
system.  Details  of  the  modifications  made  to  the  particular  data  set  used 
follow. 

2.1  DIGITAL  GRAPHIC  RECORDER  DATA 

The  process  of  digitizing  cartographic  information  has  been  used  for  various 
purposes  over  the  past  decades.  In  particular,  the  Digital  Graphic  Recorder 
(DGR)  has  been  employed  to  trace  the  lines  of  a  map  in  order  to  generate 
sequences  of  coordinate  pairs  representing  these  lines.  This  digital  data  set 
is  then  stored  on  magnetic  tape  and  is  available  in  this  form  for  computer 
processing.  Specifically,  we  are  Interested  in  the  problem  of  determining  a 
computer  internal  terrain  surface  representation  from  this  information. 

An  example  of  a  data  tape  generated  by  a  DGR  was  made  available  by  the  Defense 
Mapping  Agency  Hydrographic/Topographic  Center  (DMAHTC)  together  with  a  graphic 
plot  derived  from  this  information.  The  information  had  been  generated  from  a 
1:24,000  map  of  the  Mustang  Mountain  area  in  Fort  Huachuca,  Arizona,  —  an  area 
which  exhibits  a  large  range  of  elevations  and  slopes.  The  digitizer  resolution 
used,  0.01  inches,  corresponds  to  a  horizontal  distance  of  20  feet. 

Contour  related  data  were  extracted  from  this  cape  and  the  information  was 
reformatted  and  installed  in  a  VAX  system.  Two  basic  problems  were  encountered 
in  this  effort  concerning  the  format  of  the  data  and  the  size  of  the  files. 

First  of  all,  the  records  in  a  DGR  magnetic  tape  contain  192  CDC-1700  words. 

Each  word  contains  18  bits  of  which  only  the  high  order  16  bits  are  used.  We 
were  to  Install  this  data  set  in  a  VAX  with  a  VMS  operating  system.  The  words 
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in  this  system  are  32  blta  long  and  tha  bytes  ara  interchanged  in  each  word.  To 
solve  thia  problem  a  program  called  DGR2TAPE  waa  written  baaed  on  an  example 
provided  by  DMAHTC  (85).  Thia  program  lnterchangea  the  high  order  with  the  low 
order  bytea  in  each  word  and  than  it  regroupa  the  blta  into  16-bit  worda  by 
extracting  tha  high  order  16  blta  from  each  18-bit  group. 

The  DGR  data  tape  uaed  Included  dlgltlxed  contour*,  ridge*,  drains  and 
neatlines.  In  the  original  file  layout  (DMARTC ,  85)  that  follows,  segments  of 
such  traced  lines  are  referred  to  aa  "scan  lines": 


FIRST  RECORD:  HEADER 


CONTENTS 


Sheet  Number 


45-50 

51 

52-60 

61-68 

69-192 


Scale  (1:24000) 

Series 

tUtll. 

Registration  Marks  in  binary  (in  0.01  inches) 
Comments 


DATA  RECORDS 


WORD 


CONTENTS 


Scan  Line  Flag 

Scan  Line  Number  (I.D.) 

Tjrpe  of  Data 

Y1  (in  0.01  inches) 

Z1 

* 

* 


i 

» 

i 

\ 


* 

b 

I 

> 


I 


| 

i 


For  each  line  or  line  segment  that  la  digitised,  the  following  Information  la 
recorded:  a  "Scan  Line  Starting  Flag"  (octal  200000)  that  Indicates  a  new 
segment;  a  "Scan  Line  Number M  that  Identifies  the  segment;  a  "Data  Type  Index" 
chat  classifies  the  segments  as  contours,  neatlines,  and  ridges  or  drains;  and 
one  or  more  data  points  ar  illustrated  In  the  Data  Records  Layout  above.  Each 
record  may  contain  one  or  more  traced  lines,  or  depending  on  the  slse  of  the 
segment,  It  may  contain  only  a  portion  of  the  segment.  Mote  that  a  segment  la 
any  portion  of  a  traced  line,  covering  possibly  the  whole  line,  but  more  often 
dividing  the  traced  line  into  several  non-contlguous  segments. 

The  data  file  consisted  of  9023  blocks  of  192  words  each.  Special  constraints 
on  the  amount  of  disk  space  available  and  data  accessibility  for  the  type  of 
manipulation  Intended  guided  the  design  of  the  new  file  layout.  It  was 
necessary  to  eliminate  all  the  data  not  useful  at  this  point  of  our 
investigation.  The  first  record  or  header  record  was  eliminated  and  from  the 
data  records  only  the  contour  line  segments  were  kept.  The  resulting  data 
format  is  described  below. 


DATA  RECORDS 


WORD  CONTENTS 


1 

2 

3 

4 

5 

* 


Segment  Flag 
Segment  Number  (I.D.) 

Z-value  (Elevation  of  Contour  Line) 
XI 
Y1 
* 


*  * 

*  * 


2n+2  Xn 

2n+3  Yn 
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2.2  EDITING  THE  INPUT  DATA 


An  examination  of  the  data  showed  that  the  contour  lines  extended  beyond  the 
neatlines.  Leaving  the  contours  unedited  caused  some  undesirable  triangles  to 
appear  in  the  triangulation.  These  triangles  are  very  long  and  narrow  and  can 
be  the  source  of  large  errors  when  a  surface  la  generated.  They  are  caused  by 
the  fact  that  the  VORONOI  triangulation  produces  a  convex  hull.  A  program 
called  CLOSCONT  was  written  to  solve  this  problem  and  to  join  "neighboring" 
segments.  This  program  edits  the  contours  of  the  input  data  to  Include  only 
information  that  lies  within  a  rectangular  area  that  can  be  chosen  arbitrarily. 
The  contour  points  falling  beyond  this  boundary  are  deleted,  so  that  the 
boundary  Includes  the  last  point  of  each  contour  which  was  intersecting  it.  The 
resulting  boundary  contains  points  that  are  closer  together,  eliminating  the 
undeslred  triangles.  For  the  effort  reported  here,  the  boundaries  were  chosen 
so  that  the  number  of  points  to  be  deleted  was  kept  to  a  minimum.  CLOSCONT 
then  searches  the  digital  contour  segments  that  remain  for  adjoining  segments 
(neighbors)  and  joins  them. 

The  nature  of  the  digitised  data  is  such  that  all  points  are  in  a  grid  with  a 
point  separation  of  0.01  inches.  Thus,  neighboring  points  in  a  line  are 
separated  by  either  a  vertical  space,  a  horisontal  space,  or  a  diagonal,  with  a 
maximum  separation  between  points  of  0.01414  Inches.  An  Illustration  of  a 
digitised  line  appears  in  Figure  2.  Integer  numbers  are  used  to  Identify  each 
grid  point  although  the  distance  between  points  is  given  in  hundredths  of  an 
Inch.  Therefore,  there  is  a  finite  number  of  points  representing  a  line,  it 
was  also  observed  that  there  are  suppressed  contours  in  steep  areas  where  the 
contour  density  is  high.  In  these  steep  areas  a  small  change  in  the  horisontal 
location  of  a  point  has  greater  impact  on  its  elevation  computation  <:nd  the 
resolution  of  the  digitiser  becomes  more  important.  The  effects  caused  by  this 
fact  are  discussed  in  Section  6.1.  Recall  that  the  resolution  of  the  digitizer, 
0.01  inches,  corresponds  to  a  horizontal  ground  distance  of  20  feet. 
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Figure  2.  Portion  of  a  Digitized  Line 
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3.  THINNING  BY  TOLERANCE  BAND  ALGORITHMS 


For  the  purpose  of  reducing  or  thinning  digitized  line  data,  the  second  task  of 
this  effort,  we  selected  a  tolerance  band  method.  Thii  approach  pernits  the 
linking  of  the  sampling  process  to  the  accuracy  of  the  approximation,  selecting 
more  points  in  areas  of  high  curvature  than  in  those  of  lower  ones.  We  also 
calculate  the  contour  tangents  at  most  selected  points  except  at  the  ends  of 
segments.  A  program  called  THIN  that  performs  all  these  computations  was 
developed. 

3.1  DEFINITIONS  AND  OVERVIEW 

In  general,  the  task  of  tolerance  band  methods  is  to  select  a  suitable  sample  of 
points  from  an  initial  line  .  The  points  in  the  selected  sample  are  frequently 
called  "critical  points".  Connecting  successive  critical  points  by  straight 
lines  leads  to  a  "generalized  line"  often  desired  for  car  jphic  purposes 
(Figure  3). 

The  critical  points  are  selected  in  such  a  fashion  that  the  line  points  between 
two  successive  critical  points  are  contained  in  a  rectangular  strip,  the 
"tolerance  band",  whose  bandwidth  is  determined  from  a  given  "tolerance"  £•  The 
deviation  of  the  initial  line  from  the  general  line  can  thus  be  bounded  in  terms 
of  £. 

Several  varieties  of  tolerance  band  methods  have  been  proposed  by  Lang  (69), 
Douglas  and  Poiker  (73),  Reumann  and  Witkam  (74)  and  Williams  (78,81).  For 
details  the  reader  is  asked  to  consult  the  comprehensive  survey  by  Zoraster, 
Davis,  and  Hugus  (84).  For  this  project,  we  designed  and  used  a  one-pass 
algorithm  for  the  method  of  Reumann  and  Witkam  (74)  in  the  spirit  of  Williams 
(78).  It  will  l.  described  below. 


Figure  3.  Covering  of  an  Initial  Line  by  Tolerance  Bands  .  Solid  Dots 

Represent  Critical  Points.  Their  Heavy  Straight-Line  Connections 
Define  the  Generalized  Line. 


3.2  A.0NE-PA3S  VERSION  OF  THE  1EUMANN  AND  WITKAM  METHOD 

By 


(3.2.1)  4  1“1»***»  n»  n  >  1, 

we  denote  the  sequence  of  points  In  the  initial  line.  In  what  follows,  we  will 
describe  an  algorithm  that  implements  the  method  of  Reumann  and  Witkam  (74). 
Inherent  in  this  method  is  the  convention  that  the  first  Initial  point  P^  and 
the  last  initial  point  Pn  are  automatically  selected  as  critical  points.  In 
addition  we  impose  the  condition  that  there  should  be  no  ’’doubling  back’’  by  the 
initial  points  within  a  tolerance  band,  that  is,  '.hat 

(3.2.2)  distances  from  the  first  (critical)  point  to  subsequent  points  i 
the  tolerance  band  do  not  decrease  when  considered  in  sequence. 


Furthermore  we  require  that 

(3.2.3)  two  critical  points  have  at  least  distance  E  from  each  other, 
unless  they  represent  the  first  and  last  initial  points. 

With  these  conventions,  the  algorithm  proceeds  as  follows.  The  point  P^  is 
automatically  selected  as  the  first  critical  point.  As  to  subsequent  points  Pj., 
we  consider  the  following  conditions: 


(3.2.4)  (i) 

(ii) 
(lii) 
(iv) 


i  <  n 

distance  (Pj,  P^)  <  £ 

distance  (Pj_,  PA)  <_  distance  (Pj_,  ) 

there  exists  a  tolerance  band  of  the  second  kind  anchored 

at  P^  which  contains  the  points  P^,  1  k  i+1. 


The  point  P^  is  then  noncritical  if,  in  terms  of  the  conditions  (3.2.4), 


(3.2.5)  (i)  and  ((ii)  or  ((ill)  and  (iv))). 


The  algorithm  checks  P2>  successively  until  a  first  critical  point  Pj  is 

encountered.  If  j<n,  that  critical  point  will  play  the  role  of  P^,  that  is,  it 
will  anchor  the  search  for  a  third  critical  point,  and  so  on.  In  Figure  4  we 
give  a  pseudo-code  description  of  the  algorithm. 

The  major  part  of  the  algorithm  concerns  the  verification  of  the  tolerance  band 
condition.  Return  to  the  two  points  P^  and  ?2»  and  assume  that  n>2  and  that  ?2 
satisfies  the  C-separation  ccndition  (3.2.4.ii).  Then  there  exist  unique 
left-most  and  right-most  tolerance  bands  of  width  2  anchored  at  which  also 
contain  P2  (Figure  5).  The  flanks  of  these  tolerance  bands  are  tangent  to  the 
g-circle  around  P1.  In  particular,  the  right  flank  of  the  right-most  tolerance 
band  touches  that  circle  at  point  R,  whereas  the  left  flank  of  the  left-most 
tolerance  band  touches  at  point  L.  We  are  particularly  interested  in  the  area 
which  lies  to  the  right  of  the  line  from  R  to  L,  and  which  is  bounded  by  rays 
extending  from  the  points  R  to  L  in  the  direction  of  the  right-most  and 
left-most  flanks,  respectively  (Figure  6).  If  the  point  P3  satisfies  the 
no-doubling-back  condition  (3.2.4.iii),  then  the  subsequent  condition(3.2.4.iv) 
holds  clearly  if  and  only  if  P3  lies  in  this  area. 

In  general,  suppose  that  Pj  is  the  last  critical  point  to  have  been  determined, 
and  that  P^,  i<n,  as  well  as  the  points  between  it  and  Pj,  are  contained  in  at 
least  one  suitable  tolerance  band.  Then,  there  are  left-most  and  right-most 
tolerance  bands  —  in  extreme  cases  they  may  coincide  —  which  contain  those 
points  also.  Suppose  further  that  is  not  closer  to  Pj  than  P^  is.  Then, 

there  exists  a  tolerance  band  containing  Pj^,  and  all  previous  points  back  to 
Pj,  if  and  only  if  P^^  lies  in  the  area  indicated  in  Figure  6. 

In  this  latter  case,  there  exiBt  again  a  left-most  and  a  right-most  tolerance 

band,  each  also  covering  the  extended  set  ■(  Pj,  Pj+i .  P^,  P^+\  These 

extremal  positions  are  determined  by  the  location  of  P^+i  with  respect  to  the 
original  extreme  tolerance  bands.  To  explain  the  situation,  it  is  advantageous 
to  introduce  the  notion  of  a  "tolerance  strip"  associated  with  a  tolerance  band. 
By  this  we  mean  the  infinite  straight  continuation  of  the 
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tolerance  band  in  its  direction.  The  tolerance  strip  thus  is  bounded  by  the 
origin-end  of  the  tolerance  band  and  the  two  infinite  rays  that  extend  the 
flanks  of  the  band.  We  call  those  rays  again  the  "flanks”  of  the  tolerance 
strip.  They  are  divided  symmetrically  by  the  "center  line".  In  addition,  we 
call  the  area  described  in  figure  6  the  "area  of  flexibility"  of  with  respect 
to  the  anchor  point  Pj.  This  area  of  flexibility  turns  out  to  be  the  convex 
hull  of  the  two  tolerance  strips  associated  with  the  right-most  and  left-most 
tolerance  bands,  respectively. 

The  point  Pi+x  now  lies  in  (1)  both  tolerance  strips,  (2)  the  right-most 
tolevance  strip  only,  (3)  the  left-most  tolerance  strip  only,  or  (4)  none.  Two 
of  these  four  cases  are  illustrated  in  Figures  7  and  8.  In  case  (1),  the 
original  extreme  tolerance  strips  remain  unchanged.  However,  the  tolerance 
bands  cut  from  these  strips  are  now,  in  general,  longer  since  Pi+x  is  to  lie  on 
their  destination  ends.  In  case  (1),  only  the  right-most  tolerance  strip 
remains  the  same,  as  the  direction  of  the  left-most  strip  must  now  be  changed  to 
cover  Pi+1*  The  smallest  adjustment  that  will  achieve  this  is  a  rotation  of  the 
left-most  strip  until  its  right  flank  meets  P^+i*  This  then  characterizes  the 
new  position  of  the  left-most  tolerance  strip.  Case  (3)  is  like  case  (2), 
except  that  the  roles  of  the  right-most  and  left-most  tolerance  strips  are 
reversed.  In  case  (4)  finally,  both  tolerance  strips  need  to  be  adjusted.  The 
adjustment  is  such  that  the  right  flank  of  the  new  left-most  strip  and  the  left 
flank  of  the  new  right-most  strip  both  Intersect  at  Pj+x» 

After  a  critical  point  Pj,  J<n,  has  been  established,  and  P ^  has  at  least 
distance  £  from  Pj,  the  left-most  and  right-most  tolerance  strips  are  first, 
established  as  if  they  were  arising  from  case  (4)  above.  This  then  establishes 
the  area  of  flexibility.  Each  time  a  subsequent  point  is  found  to  lie  in  this 
area,  the  tolerance  strips  are  adjusted  and  the  area  of  flexibility  narrowed^  If 
a  subsequent  point  falls  outside  the  area  of  flexibility,  then  its  predecessor 
will  be  selected  as  a  critical  point.  This  also  will  be  the  case  if  doubling 
back  occurs,  or  if  there  are  no  points  left  (compare  the  pseudo-code  description 
in  Figure  4) . 


define  sets 

LS  -  LEFT-MOST  TOLERANCE  STRIP 
RS  -  RIGHT-MOST  TOLERANCE  STRIP 
AF  -  AREA  OF  FLEXIBILITY 
define  variables 

N  -  NUMBER  OF  INITIAL  POINTS  (INPUT) 

P(I)  -  I-TH  INITIAL  POINT  (INPUT) 

EPS  -  TOLERANCE  PARAMETER  (INPUT) 

M  -  NUMBER  OF  CRITICAL  POINTS  (OUTPUT) 
if  N<2  abort 
K:-l;  J;-l;  C(l)s«P(l) 
while  J<N  do 

LS: -EMPTY;  RS: -EMPTY;  I:«J;  CRITICAL  POINT: -FALSE 
while  not  CRITICAL  POINT  do 
I :  -1+1 

if  |P(J)-P(I)|  >EPS  then 

if  P(I)  |  LS  then  ADJUST  LS 
if  P(I)  $  RS  then  ADJUST  RS 
DETERMINE  AF 

if  I-N  or  |P(J)-P(I)|  >|P(J)-P(I+1)|  or  P(I+1)  k  AF  then 
CRITICAL  POINT: -TRUE 
K:-K+l;  J:-I;  C(K):-P(J) 

M:-K 


Figure  4.  Pseudo-Code  Description  of  the  One-Pass  Reumann  and  Witkam  Algorithm. 


CASE  (3);  The  Subsequent  Point  P^+i  Lies  in  the  Right-Most 
Tolerance  Strip  of  the  Initial  Points  {Pj,  P^,...,  P^h 
but  not  in  the  Left-Most* 


Original 

left-most 

tolerance-strip 


8.  CASE  (4):  The  Subsequent  Point  P^+i  Does  Not  Lie  in  any  of  the  Two 
Extremal  Tolerance  Strips  of  the  Initial  Points  <{Pj,  Pj+i,..., 


3.3  DRBHIIKU6  OOKTOUK.  TARGXVTS 


Consider  three  subsequent  contour  points: 

<*i»yi>»  <*2»y2>.  <*3»*3>- 

We  define  the  tangent  to  the  contour  at  the  second  point  as  the  tangent  to  the 
circle  through  these  three  points  (Figure  9).  A  vector  normal  to  the  tangent  at 
(xp.yj)  then  ij  given  by 


/ 

*12  *32 

*12  *32 

(3.3.1)  j 

;  “ 

V 

r12  r32 

r12  r32 

where 


*12  “  *1  -  *2*  *32  ’  *3  ~  *2* 

*12  "  *1  "  *2’  *32  "  *3  ’  *2* 

r12  "  *12  +  *12*  r32  *  *32  +  *32* 

This  formula  is  used  to  calculate  an  aprlorl  estimate  of  the  normal  to  the 

contour  tangent  at  each  point  that  lies  in  the  interior  of  a  contour  segment. 
Note  that  the  terrain  gradient  Is  parallel  to  this  normal.  The  same  will  hold 
for  the  gradient  of  the  subsequently  generated  synthetic  surface.  This  will 
ensure  a  more  faithful  representation  of  the  terrain. 


At  the  endpoints  of  contour  lines,  the  tangent  was  left  undetermined.  Also  if 
the  ':hree  points  through  which  the  circle  is  to  be  passed  form  an  acute  angle  at 
the  second  point,  then  the  above  determination  of  the  contour  tangent  has  little 
meaning.  In  most  of  the  later  runs,  the  tangent  was  therefore  left  undetermined 
at  such  points. 
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4.  VORONOI  -TftXAHGULATIOK  OF  URGE  SETS  OF  CONTOUR  POINTS 


After  a  suitable  sample  of  the  contour  points  hat  bean  selected,  the  nap  area  la 
partitioned  into  triangles  whose  corners  are  aanple  points.  We  speak  of  a 
" triangulation**  of  the  nap  area.  As  will  be  described  in  the  next  section,  a 
synthetic  surface  will  be  generated  by  defining  a  surface  patch  or  "element"  in 
each  triangle  of  the  triangulation.  For  this  purpose,  a  triangulation  method 
that  captures  the  proxlnity  relationships  In  the  aanple  set  is  desired.  One 
such  method  is  the  Voronoi  triangulation,  also  called  Delaunay  triangulation, 
which  is  obtained  fron  the  Voronoi  diagram  by  dualisation  in  a  fashion  described 
below. 

4.1  DEFINITIONS  AND  ALGORITHMS 

Consider  a  finite  set  S  of  points  in  the  plane.  For  any  point  P  in  S  the 
"Voronoi  polygon  of  P  relative  to  S"  is  the  set  of  all  points  in  the  plane,  such 
that  P  is  as  close  to  any  point  in  this  set  as  is  any  other  point  in  S.  The 
Voronoi  polygon  of  a  point  P  in  S  is  the  intersection  of  the  half-planes  which 
contain  P  and  which  are  determined  by  the  perpendicular  bisectors  of  the  line 
segments  connecting  P  and  the  other  points  in  S.  Thus,  the  Voronoi  polygon  of  a 
point  P  is  a  convex  polygon,  possibly  unbounded,  which  contains  P  in  its 
interior.  Given  an  arbitrary  point  X  in  the  plane  and  the  Voronoi  polygons 
associated  with  a  set  S,  then  one  and  only  one  of  the  following  statements  is 
true: 


(1)  X  lies  in  the  interior  of  one  and  only  one  Voronoi  polygon. 

(2)  X  lies  in  the  interior  of  an  edge  shared  by  two  Voronoi  polygons. 

(3)  X  is  a  vertex  of  three  or  more  Voronoi  polygons. 


It  follows  that  the  Voronoi  polygons  of  the  points  in  S  cover  the  plane  without 
overlapping,  that  is,  without  common  interior  points.  The  union  of  their  edges 
forms  a  diagram,  the  "Voronoi  diagram  for  S",  which  partitions  the  plane  into 
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tha  Voronoi  polygons.  A  point  that  lies  In  the  plane  and  satisfies  condition 
(3)  above  1«  said  to  be  a  "vertex"  of  the  Voronoi  diagram  for  S.  It  Is  called  a 
"degenerate  vertex"  whenever  It  la  a  vertex, of  more  than  three  Voronoi  polygons. 

Figure  10  Illustrates  a  Voronoi  diagram  that  hos  a  degenerate  vertex. 

The  "dual  Voronoi  diagram"  for  a  finite  set  S  of  points  in  the  plane  is  the 
diagram  obtained  by  connecting  with  straight-line  segments  those  pairs  of  points 
In  S  whose  Voronoi  polygons  relative  to  S  have  an  edge  In  common.  Figure  11 
shows  how  such  a  dual  diagram  Is  obtained  from  a  Voronoi  diagram.  The  dual 
diagram  defines  a  collection  of  non-overlapping  convex  polygons  which  cover  the 
convex  hull  of  S.  Since  each  edge  of  the  Voronoi  diagram  is  a  line  segment 
whose  end-points  are  vertices  of  the  Voronoi  diagram,  it  follows  that  there  is  a 
one-to-one  correspondence  between  the  vertices  of  the  Voronoi  diagram  and  the 
polygons  determined  by  the  dual  diagram.  In  fact,  it  follows  from  the 
definition  of  a  Voronoi  polygon  that  a  vertex  of  the  Voronoi  diagram  is 
equidistant  from  the  points  in  S  that  are  vertices  of  the  corresponding  polygon 
determined  by  the  dual  diagram,  and  it  is  closer  to  these  points  than  It  Is  to 
any  other  point  in  S.  In  general,  most  vertices  of  the  Voronoi  diagram  are 
non-degenerate,  so  that  most  of  the  polygons  defined  by  the  dual  diagram  are 
triangles.  In  the  presence  cf  degenerate  vertices,  the  corresponding  polygons 
determined  by  the  dual  diagram  can  be  partitioned  into  non-overlapping  triangles 
by  introducing  suitable  diagonals.  Thus  a  "Voronoi  triangulation"  results  from 
the  dual  diagram  for  the  set  S.  We  note  that,  while  the  dual  Voronoi  diagram  is 
uniquely  determined,  there  are  several  compatible  Voronoi  trlangulatlons  in  the 
presence  of  degeneracies.  Figures  12  and  13,  respectively,  illustrate  the  dual 
diagram  and  a  Voronoi  triangulation  that  results  from  it. 

J.  Bernal  and  S.  E.  Howe  of  the  National  Bureau  of  Standards  (NBS)  have 
generalized,  extended  and  combined  algorithms  by  Bentley,  Weide  and  Yao  (80), 
and  Bovver  (81),  to  obtain  an  algorithm  which  constructs  the  Voronoi  diagram 
and,  therefore,  a  Voronoi  triangulation,  in  "linear  expected  time"  for  a  set  S 
of  points  distributed  uniformly  in  the  interior  of  a  rectangle  in  the  plane. 

This  material  is  being  readied  for  publication  under  the  title  "Expected  0(N) 

and  0(N^^)  Algorithms  for  Constructing  Voronoi  Diagrams  in  Two  and  Three  f 
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■nw  r j  r.  r.  TW-nuuntJimntB  mi  da  kx  ru\  nn  *n  *.*,  x.n  k.iaaajtk.'wi  .w.wjviva  ’ 


Figure  12.  The  Dual  Voronoi  Diagram  or  Delaunay  Diagram  Representing  Neighbor 
Relations  Between  Points. 
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Figure  13.  Voronoi  Triangulation  Obtained  by  Adding  a  Diagonal  to 
a  Non-Triangular  Cell  in  the  Dual  Diagram  in  Figure  12. 


26 


Dimensions"  by  J.  Bernal  and  S.  E.  Howe.  A  draft  copy  is  available  upon 
request. 


In  the  resulting  algorithm,  the  set  S  which  is  to  be  triangulated  is  enclosed 
into  a  rectangle.  The  algorithm  consists  of  three  steps.  The  first  step 
divides  the  rectangle  into  approximately  N  (number  of  points  in  S)  equally  sized 
square  cells,  and  assigns  each  point  in  S  its  proper  cell.  On  the  average  there 
will  be  one  point  per  cell,  although  there  may  be  empty  cells  as  well  as  cells 
with  more  than  one  point.  Each  cell  is  of  the  form 

■{(*,  y):  x  x  <  x+b,  y  S.  y  —  y+bh 

for  some  x,  y,  b,  and  a  point  in  S  is  assigned  that  cell  if  it  lies  in  that 
area.  Cells  within  a  distance  of  two  cells  from  the  boundary  of  the  rectangle 
are  called  "outer  cells"  and  all  others,  “inner  cells." 

The  second  step  constructs  the  Voronol  polygons  of  the  points  belonging  to  inner 
cells.  Given  a  point  P  in  an  inner  cell,  a  search  for  other  points  in  S  is 
conducted  through  each  of  the  layers  of  cells  surrounding  P.  This  search 
procedure,  called  a  "spiral  search",  starts  with  the  cell  that  contains  P,  and 
then  proceeds  in  outward  direction  to  each  of  the  layers  of  cells  surrounding 
this  cell.  The  Voronol  polygon  of  P  is  progressively  built  by  intersecting  the 
half  planes  which  contain  P  and  which  are  determined  by  the  perpendicular 
bisectors  of  the  line  segments  connecting  P  and  the  points  in  S  found  through 
the  search.  A  geometrical  test  is  available,  which  permits  to  ascertain  whether 
the  Voronol  polygon  has  achieved  its  final  form.  In  most  cases  the  Voronol 
polygon  of  P  is  obtained  after  examining  only  a  small  number  of  cells  and 
points. 

The  third  step,  finally,  builds  the  Voronol  polygons  of  points  in  the  outer 
cells  by  applying  a  modified  version  of  Bowyer's  Insertion  algorithm  to  this  set 
uf  points . 
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4.2  IMPLEMENTATION 


J.  Bernal  and  S.  E.  Howe  have  Implemented  the  above  algorithm  on  a  Control  Data 
Cyber  205  at  NBS.  The  Implementation  consists  of  about  5,000  FORTRAN 
statements.  It  requires  as  input  a  tolerance  8  >  and  a  list  of  the  x  and  y 
coordinates  of  the  points  in  the  set  for  which  a  Voronoi  diagram  is  desired. 

This  list  of  points  must  be  free  of  duplication,  that  is,  the  distance  between 
any  two  points  must  always  be  above  the  tolerance  8 •  The  execution  of  the 
package  requires  approximately  34N  words  of  memory,  where  N  is  the  number  of 
points  to  be  triangulated.  Versions  of  this  package  were  successfully  trans¬ 
ferred  first  to  a  VAX  11/750  and  then  to  a  VAX  11/780  at  ETL.  This  required 
adaptations  which  are  described  below. 

One  of  the  main  objectives  of  our  work  was  to  demonstrate  the  feasibility  of 
triangulating  data  sets  of  about  40,000  to  70,000  points.  Two  difficulties 
arose  in  the  course  of  the  demonstration.  The  first  difficulty  was  that  of 
furnishing  the  package  with  the  ability  to  deal  with  data  sets  whose  members  are 
not  all  necessarily  distinct.  The  second  difficulty  had  to  do  with  memory 
restrictions  that  would  not  allow  the  execution  of  the  package  for  more  than 
50,000  points  at  a  time. 

The  "check  for  duplication”  posed  a  difficulty  because  it  was  not  possible  to 
consider  every  pair  of  points  in  view  of  their  large  number.  In  the  case  of  a 
40,000  point  set  this  would  have  amounted  to  the  examination  of  about 
800,000,000  pairs.  Fortunately,  it  was  discovered  that  a  portion  of  the  Voronoi 
triangulation  package  already  provided  a  useful  tool  for  the  solution  of  this 
problem,  namely  the  cell  structure  set  up  by  the  implementation  of  the  first 
step  of  the  algorithm.  Thu.*,  a  procedure  was  developed  which  takes  advantage  of 
this  structure  to  check  for  the  duplication  of  points.  Given  a  point  P  in  a  set 
S,  use  a  spiral  search  through  each  of  the  layers  of  cells  surrounding  P  to 
search  for  other  points  in  S.  Eliminate  any  point  found  through  the  search 
whose  distance  from  P  does  not  exceed  the  tolerance  £. 
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Terminate  the  search  as  soon  as  every  cell  that  intersects  the  closed  circular 
disk  with  center  P  and  radius  has  been  searched.  In  general,  the  disk  defined 
above  is  contained  in  the  cell  that  contains  P,  so  that  this  is  the  only  cell 
that  has  to  be  searched.  This  procedure  was  incorporated  into  the  package,  and 
was  able  to  identify  and  remove  duplications  for  data  sets  of  25,000  points  in 
less  than  30  seconds  of  GPU  time. 

The  restrictions  on  available  memory  mentioned  above  required  the  development  of 
a  "decomposition  procedure"  which  allows  the  separate  triangulation  of  a  finite 
number  of  subsets  of  a  data  set  in  such  u  way  that  the  correct  total  triangula¬ 
tion  results.  This  procedure  will  now  be  described. 

Given  a  positive  integer  k  we  telect  numbers  xq>  x^,  yo>  ?l»  8Uch  tliat  *0  ^  xk» 
y0  <  yx,  and  the  rectangle 

r  ■  *t(*»  y):  *o  £  x  £  *k»  yo  £  y  £  yi>» 

contains  the  data  set.  We  select  numbers 

*01, ’  X1L»  xl»  X1R*  x2L»  x2»  x2R»  * • • »  Xk-1,L»  xk-l»  Xk-1,R»  xkR» 
such  that 


Xq<X^<.  .  and 

X0"X0L*  x1L<x1<x1R>  x2L^x2<'x2R*  *  *  *  *  xk-l,l/xk-l<xk-l,R»  xkR*xk* 
For  each  i,  i-1,  2 . k,  we  define  rectangles  R^  and  R^  by: 

-  «K*.  y):  £ x  £  xi,r»  yo  i  y  £  yl>* 

ra  -  <(*,  y)s  x^!  <  x  <  Xi,  yo  £  y  1  yi 

k 

It  follows  that  R^  C  R^’  for  each  i,  i-1,  2,...,  k,  and  R  -  ^U^R^ 
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We  also  assume  that  the  points  (xq,  JTq).  (*o»  ^1^*  <xl»  ^0^  (xl> 

(xk-l»  y0>*  (xk-l»  yl>*  <xk*  7o)*  (xk*  yl>»  that  ia*  the  corner  Points  of  the 
rectangles  R^,  1*1,...,  k,  belong  to  the  data  set  (Figure  14).  Furthermore, 
we  define  to  be  the  set  of  points  In  the  data  set  that  belong  to  R^'  for 
each  1,  1*1,...,  k.  Assuming  that  k  and  the  numbers  xq^,  x1l,  ,  xlR,..., 

L,  xR_^,  *^-1  R,  xRR  have  been  properly  selected,  the  decomposition 
procedure  consists  of  obtaining  separate  trlangulations  T^,  k,  for  the 

sets  S^,  1*1,...,  k,  respectively  (Figure  15).  The  correct  trlangulatlon  of 

the  entire  data  set  Is  then  given  by: 


Rt  Intersects  the  Interior  of  t 


In  order  to  properly  select  k  and  the  numbers  xql,  x^,  x^,  x^r . xk-l  ,L* 

xk-1,  xR_i  R,  xRR,  a  separate  procedure  was  developed.  In  what  follows,  we 
define,  for  a  given  triangle  t  in  a  Voronoi  trlangulatlon,  x(t)  and  y(t)  to  be 
the  x,y»coordlnates  of  the  vertex  in  the  Voronoi  diagram  that  corresponds  to  t. 
Accordingly,  we  define  d(t)  to  be  the  distance  from  (x(t),  y(t))  to  any  one  of 
the  vertices  of  t.  Since  (x(t),  y(t))  Is  equidistant  from  the  vertices  of  t, 
d(t)  Is  well  defined.  In  the  following  procedure,  m  denotes  the  maximum  number 
(.  v  points  that  can  be  triangulated  with  a  single  run  of  the  package: 


Step  1.  Let  k  *  1  and  obtain  Rj,  Rj*  and  Sj.»  Let  j  *  1. 

Step  2.  If  the  number  of  points  in  Sj  does  not  exceed  m  go  to  step  3. 

Else  increase  k  to  the  next  positive  integer  for  which  *0L»  X1L*  xl» 
X1R’**’’  *k-i  l*  xk-l»  Xk-1,R'  xkR»  can  defined  with:  x0<x^<...<xR-]<xR; 
xil<x1<x1r>  *  »  xk-l ,L<xk-l<xk-l ,R*  X0*X0L »  xkR“xk»  and  the  corresponding 

Rj,  R^',  S^  ...,  k,  can  be  obtained  with  the  number  of  points  in  each 

Sj,  1*1,...,  ,  iot  exceeding  m.  Let  j  *  1. 

Step  3.  Obtain  the  Voronoi  triangulation  Tj  for  Sj.  If  j  is  equal  to  1 
let  xl*xql.  Else  define  xL  by: 

xL  *  min  -{  x(t)  -  d(t):  t  (  Tj,  2  vertices  of  t  lie  in  Rj_i> 

1  vertex  of  t  lies  in  the  interior  of  Rj  J*. 
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gulations  of  Three  Sections  Superimposed  Over  the 
gulation  of  the  Entire  Set. 


If  j  is  equal  to  k,  let  x^x^*  Klee  define  xR  by: 

xR  •  max  4  x(t)  +  d(t):  t  (  T j ,  2  vertlcea  of  t  lie  in 
1  vertex  of  t  lies  in  the  interior  of  Rj  >. 

If  Xj_j^L  does  not  exceed  x^  end  xR  does  not  exceed  XjR  go  to  step  4. 

Else, 

if  Xj_^>L  exceeds  xL  let  Xj.^^^x^,  end  if  xR  exceeds  XjR  let  XjR"xR. 

Obtain  Rj,  Rj',  Sj  end  go  to  step  2. 

Step  4.  If  j  equals  k,  stop.  Else  let  j  ■  j  +  1  and  go  to  step  3. 

This  procedure  was  Incorporated  into  the  triangulation  package  and  is  currently 
operational > 

4.3  RESULTS  OF  A  COMPUTATIONAL  EXPERIMENT 

A  Voronoi  triangulation  package  that  Includes  the  adaptations  described  above 
was  implemented  at  ETL  and  NBS.  The  feasibility  of  performing  a  triangulatlon 
with  large  sets  of  data  was  demonstrated.  For  a  specific  experiment  conducted 
on  the  VAX  11/750,  we  divided  the  triangles  into  three  classes: 

Class  1:  Those  whose  vertices  lie  exactly  on  one  contour  line  of  the  map. 
Class  2:  Those  whose  vertices  lie  exactly  on  two  contour  lines  of  the  map. 
Class  3:  And  those  whose  vertices  lie  exactly  on  three  contour  lines. 

Each  class  may  contain  triangles  that  are  intersected  by  contour  lines  other 
than  those  containing  their  vertices,  and  any  given  triangle  belongs  to  one  and 
only  one  class.* 

A  first  data  set  contained  39,645  points  and  was  decomposed  into  two  subsets.  A 
second  data  set  contained  70,249  points  and  was  decomposed  into  four  subsets. 

The  tables  below  illustrate  some  of  the  results  obtained  when  triangulating 
these  data  sets  with  the  Voronoi  triangulatlon  package.  The  times,  given  in  CPU 
seconds,  indicate  the  rate  at  which  the  package  ran  per  point. 
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TABU  rot  H  •  39,645 


Element  of 

X  Claes  1 

X  Class  2 

X  Claos  3 

Number  of 

CPU  Time 

Decomposition 

Triangles 

Triangles 

Triangles 

Triangles 

Sec/Polnt 

1 

17.20 

49.96 

32.84 

41,875 

0.044175 

2 

17.05 

48.09 

34.86 

36,451 

0.046810 

TABU  F0&  N 

-  70,249 

Element  of 

X  Class  1 

X  Class  2 

X  Class  3 

Number  of 

CPU  Time 

Decomposition 

Triangles 

Triangles 

Triangles 

Triangles 

Sec/Polnt 

1 

22.58 

59.58 

17.84 

41,134 

0.048542 

2 

22.65 

58.65 

18.70 

32,824 

0.043350 

3 

18.74 

59.53 

21.73 

33,715 

0.044793 

4 

26.00 

55.76 

18.24 

31,755 

0.052509 
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5.  SURFACE  GENERATION  WITH  CLOUGH-TOCHER  ELEMENTS 


In  this  section  we  discuss  the  construction  of  a  surface  function  passing 
through  Irregularly  spaced  given  points  for  most  of  which  both  elevations  and 
contour  tangents  are  specified.  The  planar  projections  of  all  given  points  are 
triangulated,  that  Is,  the  map  area  Is  tiled  with  triangles  whose  vertices  are 
these  points.  The  description  of  the  surface  is  in  terms  of  these  triangles:  In 
order  to  find  the  surface  elevation  for  an  arbitrary  point  in  the  map  region,  a 
triangle  containing  this  point  must  be  found.  An  evaluation  formula  or 
“element"  Is  then  invoked  using  triangle-specific  parameters.  The  "Clough-Tocher 
element"  employed  in  our  work  requires  the  determination  of  suitable  tangent 
planes  at  the  given  points.  For  this  purpose,  "local"  as  well  as  more  expensive 
"global"  methods  are  available.  A  global  method  based  on  energy  minimisation 
has  been  Implemented  and  tested  for  computational  feasibility. 

5.1  THE  CLOUGH-TOCHER  ELEMENT 

Triangulation-based  surface  interpolation  is  a  classical  computational  problem. 
Various  versions  of  the  Finite  Element  Method  (see  Zienklewicz  (71),  Blrkhoff 
and  Mansfield  (74))  are  usually  employed  in  its  solution.  The  "linear  element" 
represents  linear  interpolation  by  the  plane  through  the  vertices  of  the 
triangle  at  hand.  It  yields  a  surface  of  continuous  elevation.  However,  this 
surface  is  not  smooth  since  "creases",  that  is,  tangential  discontinuities, 
occur  along  the  boundaries  of  the  triangles. 

Nonlinear  elements  are  needed  for  smooth  surface  interpolation.  A  major 
advantage  of  smooth  surfaces  is  that  their  corresponding  surface  functions  are 
uniquely  differentiable  at  each  point  of  their  domain;  in  other  words,  there  are 
unique  gradients.  The  greater  flexibility  and  information  content  of  nonlinear 
elements  also  allows  for  a  more  precise  representation  of  the  original  data  than 
that  provided  by  linear  elements,  given  triangulations  of  comparable  densities. 
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The  "Clough-Tocher"  element  (see  Clough  and  Tocher  (65),  Lawson  (72,  76,  77))  is 
a  particularly  attractive  tool  for  smooth  surface  interpolation.  It  requires 
that  at  each  vertex  1  of  the  triangle  over  which  it  is  defined,  the  elevation 
and  the  partial  derivatives  zix  and  ziy  be  given.  The  Clough-Tocher  element 
then  is  described  by  a  function  z  -  f(x,y)  on  the  given  planar  triangle.  It 
represents  a  "surface  patch”  above  this  triangle  (Figure  17,  see  also  Figure 
16).  The  surface  patch  meets  the  prescribed  elevation  as  well  as  the  prescribed 
derivatives  at  each  vertex.  It  is  fully  defined  by  these  quantities,  in  other 
words,  by  the  three  elevations  and  the  three  tangential  planes  (gradients). 

The  following  considerations  concern  the  construction  of  an  entire  surface  from 
such  surface  patches.  In  order  to  ensure  continuity  of  elevation  between 
adjacent  triangles,  the  Clough-Tocher  element  satisfies  the  following 

(5.1.1)  Cubic  Boundary  Condition:  Along  each  triangle  edge,  the 

Clough-Tocher  element  agrees  with  a  cubic  (degree  3  or  less) 
polynomial  in  terms  of  a  variable  linearly  traversing  the  edge. 

A  cubic  polynomial  in  one  variable  is  completely  determined  by  two  elevations 
and  two  derivatives  (in  the  direction  of  the  edge).  Since  the  tangential  planes 
at  the  vertices  agree,  so  do  the  derivatives  in  the  direction  of  the  edge. 
Therefore,  the  Clough-Tocher  elements  of  two  adjacent  triangles  determine  the 
same  cubic  polynomial  on  the  edge  they  share.  It  follows  that  not  only  the 
elevations,  but  also  the  derivatives  in  the  direction  of  the  edge  agree  along 
that  common  boundary. 


In  order  to  ensure  smoothness  across  the  triangle  boundaries,  the  Clough-Tocher 
element  satisfies  the  following 


The  edge-perpendicular  derivatives  along  an  edge  are  thus  uniquely  determined  by 
their  location  on  the  edge  and  the  derivatives  at  the  vertices  at  the  ends  of 
the  edge.  Since  the  latter  agree  for  two  adjacent  triangles,  the  edge- 
perpendicular  derivatives  along  their  common  edge  agree  also.  It  was  seen  that, 
as  a  consequence  of  the  cubic  boundary  condition,  the  edge-parallel  derivatives 
agree.  Hence,  adjacent  Clough-Tocher  elements  share  tangential  planes 
everywhere  along  their  common  boundary. 

Functions  over  triangles  are  best  expressed  in  terms  of  their  "barycentrlc 
coordinates",  also  called  "triangle  coordinates."  These  are  three  real  numbers 

,  ^2 »  ^3 , 

such  that 

(5.1.3)  Xt  +  X2  +  X3  -  1 

X^x-^  +  X2x2  +  ^3X3  ■  x 

+  x2?2  +  X3y3  "  y» 

where  x,  y  are  the  planar  coordinates  of  the  point  in  question.  The  barycentrlc 
coordinates  are  functions  of  these  planar  coordinates.  To  express  these 
functional  relationships,  we  use  Zienkiewicz  notation? 

xij  "  xi  "  xj»  yij  =  yi  "  yj*  i,j  “  2’  3‘ 

We  then  have 

(5.1. A)  Xx  -  X1(x,y)  ^  (y23*x  +  x32*y  +  x2y3  “  y2x3^° 

x2  -  X2(x,y)  -  (y31*x  +  xl3*y  +  X3yx  -  y3xl 

X3  =*  X3(x,y)  -  (y12«x  +  x2,*y  +  xjy2  -  y1x2)/D, 
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where  the  denominator  is  given  by  the  determinant 


D 


1  1  1 


?1  ?2 


We  also  note  that  the  partial  derivatives  of  the  barycentric  coordinates  with 
respect  to  x,  y,  are  given  by 


(5.1.5)  ■  y23/D,  X2x  “  y3i/D,  ^3x  “  7\2^* 

xly  "  x32/D>  X2y  *  x13/°»  \)y  “  *21/D. 

The  advantage  of  barycentric  coordinates  lies  in  the  symmetric  way  the  vertices 
of  the  triangle  are  treated.  Also,  their  signs  indicate  immediately  whether  the 
point  (x,y)  lies  inside  or  outside  the  triangle:  a  negative  barycentric 
coordinate  indicates  that  the  point  lies  outside  the  triangle.  If  all 
barycentric  coordinates  are  positive,  then  the  point  lies  in  the  interior  of  the 
triangle.  On  the  boundary  at  least  one  barycentric  coordinate  vanishes. 

Vertices  are  characterized  by  single  nonzero  barycentric  coordinates  of  value  1: 


(1,  0,  0)  ,  (0,  1,  0)  ,  (0,  0,  1). 

The  barycenter  of  the  triangle  or  "centroid"  is  given  by 
(  1/3,  1/3,  1/3  ). 

It  defines  what  we  call  a  "barycentric  partition"  (Figure  18)  of  the  triangle 
(5.1.6)  T  *  B^  u  B2  u  B-j  U  Bq, 

where 
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Bj_  ■  ! 

M  ^  ^2» 

M  £  ^3^ 

Bj  m  {(X|,X2,X3). 

X2  ^  X3  | 

^2  —  ^1 ^ 

®3  “  {(Xi»X2,X3): 

X3  <  Xj, 

x3  1  ^2  ^ 

B0  ■  { "  ^2  “  ^3  *  1/3J. 

In  each  of  the  major  triangle  regions  B^,  i  ■  1*2,3,  the  corresponding 
barycentric  coordinate  is  dominated  by  the  remaining  ones: 

*  min  <(  X^ ,  X2 ,  Xj 


Each  function  representing  a  Clough-Tocher  element  is  a  cubic  polynomial  of  the 
barycentric  coordinates  in  each  of  the  major  regions  B^  of  the  barycentric 
partition  of  the  triangle.  The  function  is  continous  and  smooth  at  the 
boundaries  of  these  regions  of  the  barycenter.  We  call  such  a  function 
"piecewise  cubic"  with  respect  to  the  barycentric  partition. 

In  his  seminal  work,  Lawson  (76)  Introduces  three  "correction  functions" 

(5.1.8)  “  ^(X^  »X2 *Xj) ,  i  ■  1,2,3, 

with  which  to  describe  the  Clough-Tocher  element.  They  are  piecewise  cubic  with 
respect  to  the  barycentric  partition,  and  given  by 


9 1 


X1X2X3  +  5/6X^3  -  1/2XX2 

for 

(x^ ,x2,x3) 

-1/6X23  +  1/2X22X3 

for 

(x^,x2,x3) 

-1/6X33  +  1/2X32X2 

for 

(x^,x2,x3) 

1/81 

for 

(X|,X2,X3) 

6  Bt 
e  B2 
6  B3 

«  Bo 
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^  5/6X2  I/2X2 

for 

0^1 *^2*^3)  *  ®2 

-1/6X33  +  1/2x3^ 

for 

(X^  X2,X3)  c  B3 

-1/6XJ3  +  i/2x12x3 

for 

(H*^2»M^  €  B1 

1/81 

for 

(X^,X2,X3)  6  Bq 

x^x2x3  5/6X33  "  i/2x32 

for 

(*•1*^2  A3)  t  B3 

-1/6XJ3  +  l^X^Xj 

for 

(^>^2*^3)  ^ 

-1/6X23  +  1/2X22X1 

for 

(Xi,X2,X3)  €  B2 

1/81 

for 

Al»X2,X3)  €  Bq. 

The  motivation  for  the  choice  of  the  correction  functions  is  given  in  Lawson 
(76).  The  partial  derivatives  of  these  correction  functions  are: 


(5.1.9) 


2 

(X2X3+5/2X^  -X1)y23  +  ^1^31  ^1^2^12 

in 

B1 

(-l/2X22+X2X3)y31  +  l/2X22y12 

in 

®2 

(-l/2X32+X3X2)y^2  +  l/2X32y3^ 

in 

b3 

-l/18y23 

in 

B0 

2 

(X2X3+5/2X^  -X^)xj2  ^  X^X^x^j  ^ 

in 

B1 

(“l/2X2^X2lj)^3  l/2X2^X2i 

in 

B2 

(-l/2X32+X3X2)x21  +  l/2X32x13 

in 

®3 

D‘?2x 


D*P2y 


<X3X1+5/2X22-X2)y31  +  ^1^2^12  *  ^2^3^23 
(-l/2X32+X3X1)y12  +  l/2X32y23 

(-l/2X12+X1X3)y23  +  1 /2X!2y12 

-i/i8y3l 

(*3^l'*‘5/2X2  “^2^x13  ^1^2*21  ^  ^2^3X32 

(-1/2X32+X3X1)x21  +  1/2X32x32 

(-1/2X12+X1X3)x32  +  1/2X12x21 
1-1/18x13 


D'P3x 


'(XiX2+5/2X32-X3)yi2  +  ^2^3^23  +  ^3xly31 
(-l/2X12+X1X2)y23  +  l/2X12y31 
(-l/2X22+X2X1)y3l  +  l/2X22y23 

-1/l8y12 


2 

(^1^2+^^3  — ^2^3X32  ^3^1^13 

(-l'~  2+X1X2)x32  +  1/2X12x13 

(-l/2X22-fX2X1)x13  +  1/2X22x32 


1-1/18: 


In  B2 
In  B3 
In  B^ 
In  Bq 
in  B2 
in  B3 

in  Bt 
in  Bq 

in  B3 
in  Bl 
In  B2 
in  B0 
in  B3 
in  Bx 
in  B2 
in  Bq. 
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Now  let 


(5.1.10)  Mjt  -  *ixXji  +  *iyyji»  A»j  "  *1  2»  3» 

These  six  quantities  represent  the  directional  derivatives  at  vertex  1  with 
respect  to  the  direction  vector  (Xj^.yjj),  which  represents  a  directed  edge  of 
the  triangle.  It  will  be  convenient  to  use  the  additional  abreviations 

(5.1.11)  “  l/2(Mj^  +  M^j),  “  l/2(Mjj  ~  M^j)  _  z^. 

Note  that  for  a  linear  function  z  -  f(x,y),  -  -M^j  ■  Zj^,  and  for  a 

quadratic  function,  Mjj  -  -  ^zji*  As  a  consequence,  the  coefficients 

and  vanish  for  linear  functions,  and  the  coefficients  C^j  for  quadratic 
functions.  We  also  denote  by 

(5.1.12)  Lt,  i  -  1,  2,  3, 

the  euclidean  length  of  the  edge  opposite  to  vertex  i.  The  Clough-Tocher 
element  is  now  of  the  form: 

(5.1.13)  z  -  \lzl  +  X2z2  +  X3Z3 

+  Q23^2x3  +  ^31X3X1  +  ^12X1X2 
+  C23V1  +  C31V2  +  C12V3, 

where  the  functions  V^,  i  -  1,2,3,  are  given  by 

“  X2X3(X2~X3)  +  13(L2+L3)(L2~L3)/Lj^]p^  "  p2  +  ^3 

V2  “  X3X^ (X3~X^ )  +  [3(L3+L^)(L3_L^)/L2^lp2  -  P3  + 

V3  “  X^X2(X^-X2)  +  [3(L^+L2)(L^~ L2)/l>3^]p3  ~  +  p2* 

Note  that  the  expression 

X^z^  +  X2z2  +  X3Z3 
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represents  the  linear  element,  namely  the  plane  passing  through  the  elevations 
given  at  the  vertices  of  the  triangle.  The  portion 


Q23X2X3  +  Q31.X3X1  +  Ql2xlx2 


is  a  quadratic  function  in  the  entire  triangle.  Note  that  quadratic  functions 
satisfy  the  Linear  Derivative  Condition  (5.1.?.)*  The  remaining  portion  is  a 
piecewise  cubic  function  with  respect  to  the  barycentrlc  partition.  This 
function  also  satisfies  the  Linear  Derivative  Condition  as  a  result  of  the 
choice  of  correction  functions  p^. 

It  is  easy  to  derive  various  expressions  for  the  gradient  of  the  Clough-Tocher 
function.  The  gradient  components  are  the  partial  derivatives  with  respect  to 
x,  y,  and  may  be  written  as 

(5.1.14)  D-sx  -  *iy23  +  z2y31  +  z3y12 

+  Q23(X3y31+X2y12)  +  Q3l(xly12+X3y23)  +  ^12^X2y23+Xly31^ 
+  C23D*Vlx  +  C31°'V2x  +  C12D-V3x 

D-Zy  -  *1X32  +  z2x13  +  z3x21 

+  Q23(x3x13+x2x21^  +  Q3l(xlx21+X3x32^  +  Ql2^X2x32+Xlx13^ 
+  C23D.Vly  +  C31D»V2y  +  C12D.V3y, 

where 

D*V^X  ■  X3(2X2~X3)y3^  -  X2(2X3-X2)yj2 

+  [3(L2+L3)(L2-L3)/L12]D.plx  -  D«p2x  +  ^*P3x 

D‘V2x  ”  xl(^x3“xl)y12  -  X3(2X^-X3)y23 

+  [3(L3+L^)(L3~L^)/L22]D»p2x  "  ®’P3x  +  ®*Plx 

^  ^3*  X2^xl"X2)y23  -  X^(2X2“X^)y3^ 

+  [SCLj+LjXLj-LjJ/L^JD'p^  -  D,Pix  +  d'P2x 


D*Viy  ■  X3(2X2-X3)xi3  “  k2(2Xj-X2)x2i 

+  [3(L2+L3)(L2-L3)/L12lD^ly  -  D-p2y  +  D*p3y 

D *V2y  *  X^(2X3~X^)x2^  ■  X3(2X^-X3)x32 

+  [3(L3-H.1)(L3-L1)/L22]D*p2y  "  D*?3y  +  D‘?ly 

D»V3y  ■  X2(2X^~X2)x32  -  X^(2X2~X^)x^3 

+  l3(L1+L2)(L1-L2)/L32]D*93y  -  D‘9iy  +  D‘?2y‘ 

A  more  transparent  formula  for  the  gradient  can  be  derived.  For  X^  ■  1, 
i-1,2,3,  that  is,  at  the  three  corners  of  the  triangle,  pix  -  ply  *0.  It  thus 
follows  from  (5.1.14)  that 

D#2lx  *  *1*23  +  z2y31  +  *3y12  +  ^12+C12)y31  +  ^3l“C31^y12 
d,*2x  “  zly23  +  z2y31  +  z3y12  +  ^23+C23)y12  +  (Q12_C12)y23 
D  *z3x  "  ziy23  +  z2y3l  +  z3y12  +  ^l^Sl^S  +  <Q23_C23)y31 

D*zly  *  zlx32  +  Z2X13  +  Z3X21  +  ^12+C12)x13  +  ^3l“C31)x21 
D  *z2y  “  zlx32  +  Z2X13  +  Z3X21  +  ^23+C23)x2l  +  (Q12"C12)x32 

D,z3y  *  zlx32  +  Z2X13  +  Z3X21  +  <Q31+C31)x32  +  <Q23“C23)x13» 

where  stix,  zly,  i»l,2,3,  are  the  prescribed  gradient  components  at  the  vertices 
of  the  triangle.  It  follows  that 

XlD-*ix+X2D**2x^3D*z3x  "  zly23+z2y31+z3y12 

+  Q23(X3y31+X2y12)-H)31(X1y12+X3y23)+Q12(X2y23+X1y31) 

-  C23(X3y31-X2y12)-C31(X1yl2-X3y23)-C12^X2y23-X1y31) 


X1D**ly+X2D-*2y+k3D**3y  “  *1x32+*2x13+*3*21 

+  Q23^3x1  3+^2x2l)+Q31^Xlx21+X3x32  ^"^12^X2x32"*"Xlx13 
-  C23(X3X^3“k2x2^)-C3^(k^x2^“X3X32)-Cj2(X2X32~X^x^3) « 

Formula  (5.1.14)  can  now  be  rewritten  aa 

(5.1.15)  *x  -  X1*lx  +  *2*2x  +  X3*3x  +  C23Alx  +  C3lA2x  +  C12A3x 

*y  “  Xl*ly  +  X2*2y  +  X3*3y  +  C23Aly  +  C31A2y  +  C12A3y* 

where 

Alx  ■  Vlx  +  <x3y3l-X2yi2)/D*  Aly  ”  vly  +  (X3X13“X2*21)/D 

^x  "  V2x  +  (xiyi2‘X3y23)/D»  A2y  "  V2y  +  ^Xlx21”X3x32^° 

A3x  "  V3x  +  C^2y23“Xly31)/D*  A3y  “  v3y  +  <X2X32“X1X13)/D* 

To  sum  up,  given  three  points  with  their  elevations  and  gradients,  that  Is, 
given  15  quantities 

xi»  yi»  *i»  *ix*  *iy»  1  "  1>  2>  3» 

the  Clough-Tocher  element  defines  a  surface  that  assumes  those  prescribed 
elevations  and  derivatives  (gradients).  In  order  to  calculate  the  elevation  at 
an  arbitrary  specified  point,  we 

o  compute  the  auxiliary  quantities  and  from  the  above 

15  quantities 


o  evaluate  the  correction  functions  for  the  barycentric  coordinates 
X^,  X2,  X3  with  respect  to  the  three  given  points  in  the  plane 

o  determine  the  values  and  enter  them  into  the  Clough-Tocher  formula 
(5.1.13).  The  gradient  components  are  computed  analogously  using 
formula  (5.1.15). 
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5.2  SURFACE  GENERATION  OVER  A  TRIANGULATED  REGION 


In  order  to  pass  a  smooth  surface  through  given  elevations  at  irregularly  spaced 
points,  one  may  partition  the  region  in  which  the  surface  is  to  be  defined  into 
triangles,  specify  gradients  along  with  the  evevations  at  the  vertices  of  these 
triangles,  and  generate  the  resulting  Clough-Tocher  patch  in  every  triangle.  The 
Clough-Tocher  patches  are  designed  in  such  a  fashion  that  they  lit  together 
smoothly  along  common  boundary  edges.  However,  there  is  still  a  missing  link. 

In  order  to  fully  define  our  synthetic  surface,  we  have  yet  to  provide  the 
gradients  (tangential  planes)  at  the  vertices  of  the  triangulation.  There  are 
several  approaches  to  finding  suitable  gradient  values.  One  is  to  examine 
neighboring  elevations  and  to  estimate  a  suitable  position  of  the  tangential 
plane  at  the  point  in  question  by  using  local  Interpolation  or  least  squares 
regression.  The  success  of  spline  techniques  suggests  a  different  approach.  A 
spline  can  be  interpreted  as  an  idealized  mechanical  structure  consisting  of 
"thin  beams"  which,  when  forced  through  specified  points,  assume  a  position  In 
which  a  surrogate  elastic  energy  Is  minimized.  Since  oscillatory  behavior  is 
associated  with  high  elasti ~  energy,  minimizing  elastic  energy  tends  to  minimize 
oscillations.  In  this  work,  we  extend  this  approach  to  two  dimensions.  We 
consider  a  mechnical  structure  consisting  of  thin  beams  of  equal  "thickness" 
along  the  edges.  They  are  joined  together  at  vertices  by  small  "thin  plates” 
which  represent  tangential  planes  forced  to  be  met  by  the  adjacent  thin  beams. 
The  sum  of  the  surrogate  energies  of  all  thin  beams  is  then  minimized.  This  is 
achieved  by  varying  the  positions  of  the  thin  plates  to  which  the  adjacent  thin 
beams  must  be  tangential  while  passing  through  prescribed  elevations  and 
satisfying  other  side  conditions  that  may  have  been  specified  for  selected 
vertices.  Once  the  idealized  mechanical  structure  has  found  its  optimal,  that 
is,  energy-minimal  position,  the  gradient  at  each  vertex  is  defined  by  the  tilt 
of  its  thin  plate.  The  resulting  surface  is  supported  by  thin  beams  much  as  the 
fabric  of  an  umbrella  is  spanned  by  its  ribs. 

Mathematically,  the  surrogate  elastic  energy  is  a  positive  definite  quadratic 
form  in  those  parameters  that  are  permitted  to  vary.  Setting  the  partial 
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derivatives  of  the  energy  with  respect  to  these  variables  to  sero  yields  the 
optimality  conditions  in  the  form  of  a  large  system  of  equations. 

The  classic  Gauss-Seidel  method  for  solving  a  system  of  linear  equations  is  an 
iterative  procedure  where  each  iteration  consists  of  passing  in  sequence  through 
all  the  equations  of  the  system.  Each  equation  is  used  for  determining  a  new 
value  for  a  particular  unknown  variable  while  keeping  the  other  unknowns  fixed. 
The  method  thus  requires  an  initial  value  for  each  unknown.  To  start,  the  first 
equation  is  transformed  into  a  linear  equation  of  only  the  first  unknown  by 
substituting  into  this  equation  the  initial  values  of  all  remaining  variables. 
This  equation  then  yields  an  improved  value  for  the  first  unknown  variable.  This 
value,  along  with  initial  values  for  the  third  and  subsequent  variables,  enters 
the  second  equation,  which  then  yields  an  equation  for  the  second  variable 
alone,  and  so  on. 

We  modify  che  Gauss-Seidel  procedure  slightly.  To  this  end  we  observe  that  each 
variable  of  the  system  of  linear  equations  is  "associated"  with  a  particular 
vertex.  The  associated  variables  at  each  vertex  satisfy  "local  optimality 
conditions".  These  optimality  conditions  have  a  structural  Interpretation:  They 
express  the  conditions  for  the  structural  parameters  represented  by  the 
associated  variables  to  assume  minimum  energy  values,  supposing  that  the 
structural  parameters  at  all  other  vertices  remain  fixed.  The  local  optimality 
conditions  again  take  the  form  of  linear  equations  in  the  variables  associated 
with  the  vertex.  In  a  typical  case,  the  thin  plate  at  the  vertex  is  tilted  into 
the  best  position  it  can  assume,  given  the  tilts  and  elevations  at  neighboring 
vertices.  The  local  optimality  conditions  then  define  this  locally  optimal 
tilt.  It  can  be  shown  that: 

(5.2.1)  THEOREM:  The  linear  system  of  equations  for  optimizing  the 

position  of  the  idealized  mechanical  structure  is  equivalent  to 
the  combination  of  all  local  optimality  conditions. 

Our  variation  of  the  Gauss-Seidel  method  now  is  to  pass  through  all  vertices  in 
sequence,  solving  at  each  vertex  the  local  optimality  conditions.  It  was  a 
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major  concern  at  the  start  of  this  project  whether  this  procedure  was 
computationally  feasible.  We  found  that  the  time  needed  for  a  Gauss-Seidel 
iteration  of  the  kind  described  above  was  well  within  an  acceptable  time  frame. 
As  discussed  in  Section  6,  this  marks  a  main  accomplishment  of  this  feasibility 
study . 

In  what  follows,  we  will  list  the  local  optimality  conditions  for  various  sets 
of  associated  variables.  We  distinguish  several  "types  of  vertices"  according 
to  their  kinds  of  associated  variables .  A  particularly  important  case  is  the 
one  in  which  the  tangent  to  the  contour  curve  is  given  along  with  its  elevation. 
In  this  case,  the  direction  of  the  gradient  is  given  —  it  is  perpendicular  to 
the  contour  tangent  —  and  the  only  variable  to  be  determined  is  the  length  or, 
rather,  a  positive  or  negative  "gradient  multiplier"  It  represents  the  only 
variable  associated  with  a  vertex  of  this  type. 

it 

We  use  a  three  letter  code  to  characterize  vertex  types.  The  first  letter  of 
the  code  refers  to  elevation:  it  is  'E'  if  the  elevation  is  given,  and  'N' 
otherwise.  In  the  latter  case,  the  elevation  is  an  associated  variable.  The 
second  and  third  letters  refer  to  gradient  components  z%  and  zy,  respectively. 
Letter  'X'  in  the  second  position  indicates  that  the  x-component  of  the  gradient 
is  given,  and  'N'  in  this  position  indicates  that  it  is  not.  Analogously,  letter 
'Y*  in  the  third  position  indicates  that  the  y-component  of  the  gradient  is 
given,  and  'N'  indicates  that  it  is  not.  Finally,  the  combination  ’RN’  is  found 
in  positions  two  and  three,  if  the  gradient  direction,  but  not  the  gradient 
itself,  is  specified.  For  example,  'ERN'  signals  the  case  in  which  both 
elevation  and  gradient  direction  are  prescribed,  leaving  the  gradient  multiplier 
as  the  only  associated  variable.  Any  occurrence  of  ’N’  indicates  an  unknown 
associated  variable.  In  particular  ’NNN'  is  used  if  all  three  quantities, 
elevation  and  gradient  components,  are  to  be  determined. 

The  surrogate  energy  of  a  thin  beam  of  length  L  is  given  by 

s“L 

(5.2.2)  E  -  J  z"(s)2ds . 

B"0 
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It  is  a  well-known  result  of  the  theory  of  thin  beams  that  the  above  expression 
is  minimized  by  cubic  functions  of  the  distance  s  along  the  projection  of  the 
beams  into  the  plane.  These  cubic  functions  are  uniquely  determined  by  the 
elevations  and  slopes  at  the  end  points  ("Hermite  interpolation”).  The 
following  local  optimality  conditions  are  derived  using  these  facts. 

In  describing  local  optimality  conditions  at  vertex  i,  we  let  the  vertices  j  run 
through  the  "star”  of  vertex  i,  namely  the  following  set  of  vertices  j: 

(5.2.3)  S(i)  ■  {  j  :  j  is  connected  by  an  edge  to  vertex  i>. 

The  abbreviations  (5.1.10)  and 

(5.2.4)  =* 

are  used  to  denote  the  edge-directional  derivatives  and  the  distances  between 
vertices,  respectively.  In  addition,  we  will  need  the  quantities 

(5.2.5)  ji^  -  plxxji  +  piyyji> 

where 

Fix’  Fiy«  Fix2  +  Fiy2  *  1 

are  the  components  of  the  gradient  direction  at  vert.ex  i  normalized  to  length  1. 
The  resulting  optimality  conditions  for  all  types  of  vertices,  except  type 
'EXY',  are  displayed  in  Figure  19.  'EXY'  represents  the  fully  specified  case  in 
which  there  are  no  associated  variables  to  be  determined. 

To  sum  up,  the  vertices  of  the  triangulation  are  divided  prior  to  surface 
generation  into  types  depending  on  which  surface  parameters  are  given  and  which 
have  to  be  determined.  These  types  are  described  by  the  letter  codes 

(5.2.6)  'ERN',  'ENN',  ' EXN ' ,  'ENY',  'EXY', 

'NRN' ,  'NNN',  'NXN' ,  'NNY' ,  ' NXY ' . 
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For  each  of  these  ten  types  —  except  the  fully  specified  case  'EXY'  —  there  is 

an  update  formula  by  which  the  surface  parameters  at  the  triangulation  vertices 
are  calculated  from  the  current  values  of  the  surface  parameters  at  neighboring 
vertices.  Sequentially  updating  every  vertex  in  this  fashion  constitutes  a 
Gauss-Seidel  iteration.  Such  iterations  are  repeated  until  the  changes  in  the 
variables  remain  within  a  given  tolerance  or  a  specified  limit  on  the  number  of 
iterations  is  reached. 

5.3  FINDING  THE  RIGHT  TRIANGLE 

Suppose  a  surface  is  specified  in  terms  of  triangular  Clough-Tocher  patches.  In 
order  to  evaluate  the  elevation  at  a  given  point  in  that  surface,  a  triangle 
containing  that  point  needs  to  be  found.  In  what  follows  we  describe  a  method 
for  finding  such  a  triangle.  This  method  is  intended  for  applications  in  which 
not  just  one  point,  but  sequences  of  such  points  are  given,  most  of  which  are 
moreover  close  to  each  other.  The  sets  of  sequential  contour  points  as  they 
arise  from  digitized  contour  information  are  a  case  in  point.  A  closely  spaced 
regular  grid  is  another.  In  these  cases,  a  given  point  in  the  sequence  will  lie 
with  high  probability  in  the  same  triangle  or  in  a  triangle  directly  adjacent  to 
the  triangle  of  the  previous  point. 

The  method  we  use  for  finding  a  triangle  containing  a  given  point  relies  for  its 
efficiency  on  the  above  observation.  Before  searching  for  the  triangle  of  a 
given  point,  the  method  requires  that  an  arbitrary  starting  point  be  specified 
for  which  a  triangle  containing  it  is  known.  We  then  move  from  this  starting 
point  straight  towards  the  given  point  until  we  reach  the  boundary  of  the 
starting  triangle  or  the  given  point  itself,  whichever  happens  first.  If  the 
boundary  is  reached  at  a  non-vertex  point,  that  is,  somewhere  in  the  interior  of 
a  boundary  edge,  then  the  unique  adjacent  triangle  can  be  readily  identified 
using  our  triangulation  data  structure,  and  we  continue  moving  in  that  triangle 
as  far  as  possible  or  necessary.  In  the  unlikely  case  that  a  vertex  is 
encountered,  all  triangles  adjacent  to  this  vertex  are  examined  in  sequence 
until  one  is  found  in  which  progress  can  be  made  towards  the  given  point.  The 
processes  are  repeated  until  the  given  point  is  reached.  The  last  triangle  in 
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Cage  ERN:  elevation  z,  and  gradient  direction  (MioMtv)  specified,  gradient 
multiplier  %  to  be  determined. 


x  ?  (Z)s  l3'”  + 


Cage  EXN:  elevation  z,  and  gradient  component  z,s  specified,  gradient  com¬ 
ponent  XiV  to  be  determined. 

x  2$^  1 3  =  ^  ~  +  Wi/] 


Case  ENY:  elevation  ti  and  gradient  component  ^  specified,  gradient  com¬ 
ponent  tiz  to  be  determined. 


*•*  x  ^51  (Lji)9  ~  ?  (Z^yst3^*  “  M >'] 


Case  ENN:  elevation  *  specified,  gradient  (z,*,z,v)  to  be  determined. 


(*#) 


«*xS?fi^  +  «‘x2£^  " 

^x2?^  +  *'x2?f@“?^+^> 


Case  NP.N:  gradient  direction  specified,  elevation  z,-  and  gradient 

multiplier  %  to  be  determined. 


* x  (i,0s  +  ',,x3?(^F  =  ?  (Z^1®*' +  3JWi'1 


«x3?#  +  -*x2?{z#  =  ?^+,w 


(t*ji 


Figure  19.  Optimality  Conditions  for  yertices  by  Type. 


54 


Case  NXN:  gradient  component  z^  specified,  elevation  z,  and  gradient  com¬ 
ponent  Zit  to  be  determined. 

+  '■> x  3  ^  (L^y  = 

*x,?(Mi  +  =  ?(^)*|s%  ~2***‘ +w<1 


Case  NNY:  gradient  component  tiV  specified,  elevation  a,  and  gradient  com* 
ponent  z»  to  be  determined. 

+  *“xS?(^i  =  ^  (£*)• '6*'  ~  ,v,iXi' +  3Mi/i 
*xs?(^ +  -x2?{S  - 


Case  NNN:  no  parameters  specified,  elevation  *,  and  gradient  (*«,*»)  to 
be  determined. 

l‘x6?®5+,i*x3?(i55)+I"xS^(^5  = 

-xs?(5+-x2?!S^x2?P  - 
-if,?i5+-x2?^x2?gf  - 


Case  NXY:  gradient  (z,,,  z»>)  specified,  elevation  z,  to  be  determined. 


Figure  19.  (Continued) 


the  chain  of  triangles  obtained  in  the  course  of  this  procedure  will  contain  the 
given  point  (Figure  20) . 

It  is  clear  that  this  procedure  will  work  best  for  a  sequence  of  points  in  which 
the  previous  point  can  serve  as  a  close  starting  point  for  the  task  of  locating 
the  given  point  in  a  triangle.  For  points  in  a  regular  grid  arranged  by 
sequential  rows  the  following  procedure  is  used.  Locate  the  first  point  of  the 
first  row  from  an  arbitrary  starting  point.  Make  a  note  of  the  first  row  and 
its  triangle  for  further  reference.  Use  it  also  as  starting  point  for  the 
second  point  in  the  first  row.  Then  use  the  second  point  and  its  triangle  as 
starters  for  locating  the  third  point,  and  so  on,  until  the  end  of  the  row  is 
reached.  Then  retrieve  the  first  point  and  its  triangle  and  use  them  as 
starters  for  locating  the  first,  point  in  the  second  row.  This  point  and  its 
triangle  are  again  kept  for  further  reference,  while  the  second  row  is  traversed 
in  the  same  manner  as  the  first  row.  This  process  is  repeated  for  the  remaining 
rows  • 
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6.  RESULTS  AND  CONCLUSIONS 


In  order  to  test  the  computational  feasibility  of  the  approach  proposed  in  the 
previous  sections,  a  pilot  implementation  was  applied  to  the  Mustang  Mountain 
data  set  (see  Section  2).  CPU  times  were  recorded  and  experience  about  memory 
requirements  was  gained.  As  a  preliminary  testing  procedure  the  residuals 
obtained  when  trying  to  recover  the  contour  elevations  were  gathered  and 
analyzed.  After  describing  the  set-up  of  the  experiment,  we  report  its  results 
and  the  observed  computational  effort. 

6,1  SETTING  UP  THE  EXPERIMENT 

The  pilot  implementation  consists  of  several  independent  modules  whose  output 
files  serve  as  input  files  to  subsequent  modules.  The  Interrelationship  of 
these  modules  and  their  interconnecting  files  is  schematically  described  in 
Figure  21. 

The  original  input  file  has  the  format  (see  Section  2)  of  Digital  Graphic 
Recorder  Data  (DGR)  and  refers  in  our  case  to  the  Mustang  Mountain,  Fort 
lluachuca,  area.  These  data  are  input  into  the  module  EDIT,  which  extracts  and 
edits  contour  information  as  described  in  Section  2  of  this  report.  The 
resulting  edited  digital  contour  file  is  then  sampled  and  contour  tangents  are 
determined  in  module  THIN.  The  sampling  method  is  described  in  Section  3  of 
this  report.  The  resulting  sample  is  first  fed  to  the  module  VORONOI  (see 
Section  4),  which  determines  a  Voronoi  triangulation  of  the  sample  points. 
Accordingly,  the  main  output  o*  this  module  is  a  "triangle  table"  that  lists  the 
vertices  and  neighbors  of  each  triangle.  In  addition,  duplicate  sample  points 
are  Identified  and  points  that  need  to  be  added,  such  as  map  corners,  are 
recorded.  This  information  is  utilized  by  the  UPDATE  module  which  creates  the 
final  "data  base".  This  data  base  is  needed  along  with  the  triangle  table  to 
generate  the  surface  in  module  SURFACE.  This  module  also  provides  options  for 
evaluating  the  elevation  of  either  random  points  or  points  in  a  regular 
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Figure  21.  System  Layout. 


rectangular  grid.  The  design  of  the  testing  procedure  is  based  on  the  following 
observation.  When  the  sample  points  enter  the  trlangulatlon  process  they  are 
considered  as  points  in  a  plane.  No  longer  do  they  hold  any  direct  relationship 
to  other  points  in  the  contour  as  such.  Therefore ,  these  original  contour 
points  are  in  essence  independent  of  the  sampled  ones.  The  testing  procedure 
consists  of  evaluating  the  surface  at  the  original  contour  points  and  comparing 
the  results  to  the  given  elevations,  yielding  a  set  of  residuals.  These 
residuals  provide  an  indication  of  the  ability  of  our  algorithm  to  "recover"  the 
original  data. 

The  full  original  data  set  contains  features  such  as  lake  shores,  lake 
hatchings,  dams,  and  peak  elevations,  which  the  algorithm  in  its  present  state 
of  development  does  not  yet  handle  in  an  accurate  fashion.  Indeed,  the  primary 
goal  of  the  effort  reported  here  is  to  establish  computational  feasibility.  Also 
some  limitations  of  the  testing  procedure  itself  need  to  be  pointed  out. 

First,  the  digitized  data  themselves  carry  a  "digitization  error"  so  that  it  is 
not  necessarily  desirable  to  reproduce  the  digitized  data  precisely.  Indeed,  if 
the  recovered  contours  represent  a  "smoothing"  of  the  digitized  lines,  they  may 
be  more  representative  of  the  true  surface  than  the  given  digitized  contour 
points. 

Second,  in  flat  terrain,  the  nonsampled  digitized  contour  points  tend  to  be 
close  to  the  boundary  of  triangles,  so  that  the  behavior  of  the  surface  in  their 
interior  is  monitored  to  a  less  extent  than  in  steep  terrain,  where  more  contour 
lines  cut  through  the  Interior  of  triangles.  In  order  to  retain  interior 
monitoring,  the  sample  was  deliberately  chosen  somewhat  smaller  than  indicated 
for  the  purpose  of  improved  accuracy,  but  still  large  enough  to  test 
computational  feasibility. 

Third,  the  vertical  deviations  measured  by  the  residuals  tend  to  be  dispropor- 
tionally  large  in  steep  terrain.  Ideally,  the  3-dimensional  distances  from  the 
true  data  points  to  the  generated  surface  should  be  evaluated  in  order  to 
determine  the  accuracy  of  the  latter.  In  flat  terrain  the  vertical 
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deviation,  that  la,  the  realdual,  provides  a  good  approxlaatlon  to  the  true 
distance.  However,  this  Is  not  the  case  for  steep  terrain.  For  this  reason  we 
introduce  an 

residual 

(6.1*1)  "adjusted  residual"  ■ 

\|  „2  + 1 

based  on  a  measure  m  of  "steepness"  of  terrain.  The  length  of  the  gradient 
(5.1.16)  appears  to  be  a  natural  measure  of  steepness.  However,  the  example  of 
a  mountain  peak,  where  the  slope  is  by  definition  sero,  shows  that  the  slope  at 
a  single  point  Is  not  a  good  indicator  of  steepness.  For  the  purposes  of  this 
report,  we  determine  the  triangle  containing  the  point  in  question  and  then 
chose  the  vertex  slope  of  largest  magnitude,  taking  into  account  that  the  unit 
length  In  the  plane  is  20  feet  (Section  2.1).  The  adjusted  residual  would 
represent  the  3-dimensional  surface  precisely  if  the  surface  were  linear,  that 
Is,  a  tilted  plane  In  a  suitable  neighborhood  of  the  data  point.  In  general, 
however.  It  is  still  an  approximation,  but  a  better  one  than  the  unadjusted 
residual.  For  horizontal  terrain,  m  »  0  and  the  adjusted  residual  equals  the 
original  one.  In  all  other  cases,  the  adjusted  residual  is  smaller.  In  our 
experiment,  we  collected  statistics  on  both  types  of  residuals. 

In  Mandel,  Witzgall,  and  Bernal  (86),  the  results  of  a  first  test  run  were 
reported.  For  this  run,  the  unadjusted  residuals  for  the  full  set  of  digitized 
contour  points  were  collected,  including  lake  hatchings  and  dams,  even  though 
the  algorithm  is  not  yet  equipped  to  handle  nonsmooth  terrain,  as  pointed  out 
above.  Nevertheless,  95%  of  the  residuals  were  between  +12.5  feet,  indicating 
t  \at  at  least  90%  did  not  deviate  more  than  halfway  to  the  next  contour  line . 
Furthermore,  an  analysis  of  the  biggest  residuals  led  to  the  discovery  of 
several  contour  lines  whose  altitudes  had  been  apparently  miscoded.  For  the 
purpose  of  the  more  extensive  experiments  reported  here,  the  original  data  set 
was  purged  of  lake  hatchings  and  the  altitudes  were  recoded  for  the  above 
contour  lines.  In  addition,  one  contour  line  representing  a  dam  was  removed 
from  consideration.  In  what  follows,  the  results  are  based  on  this  "sanitized" 
data  set. 
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6.2  RESULTS 


The  experiment  reported  here  comprises  three  runs,  all  for  the  same  "sanitised" 
Mustang  mountain  data  set  of  about  38,000  sample  points  giving  rise  to  about 
75,000  triangles.  For  the  first  and  main  run,  a  histogram  of  the  (unadjusted) 
residuals  is  displayed  in  Figure  22.  In  addition,  the  standard  statistical 
quantities  such  as  expected  value  ("average),  standard  deviation,  maximum  and 
average  absolute  value  are  reported,  the  latter  two  also  for  the  adjusted 
residuals  (6.1.1).  The  second  run  did  not  utilise  the  contour  tangent 
information  and  the  third  run  used  linear  rather  than  Clough-Tocher 
interpolation  on  the  given  set  of  triangles. 

The  timing  of  the  procedure  is  broken  down  by  major  steps,  including  the 
generation  of  a  901  X  901  rectangular  grid,  which  we  expect  to  be  of  a  sise 
relevant  to  prospective  applications.  The  calculations  were  timed  and  carried 
out  on  a  VAX  11/780  system  at  the  Engineer  Topographic  Laboratories.  The 
observed  CPU  times  for  the  different  steps  follow 


o  Step  1:  Editing  Digitised  Input  Data .  25  min 

o  Step  2:  Thlning  and  Tangent  Determination  ....  9  min 

o  Step  3:  Voronoi  Triangulation  and  Update  .  21  min 

o  Step  A:  Surface  Generation .  20  min 

o  Step  5:  Grid  Determination  (901  X  901)  .  26  min 

Total  CPU  Time  . .  101  min 


The  timing  of  Step  5  represents  an  improvement  over  the  one  reported  in  Mandel, 

Uitzgall,  and  Bernal  (86).  i 

As  discussed  before,  the  testing  procedure  consisted  of  calculating  and 

analysing  the  residuals  of  the  elevation  of  the  original  contour  points.  The  j 

histogram  of  these  residuals  is  displayed  in  Figure  22.  Other  statistics  { 

calculated  include:  I 


FREQUENCY  (K) 


HISTOGRAM  OF  RESIDUALS 


Figure  22.  Histogram  of  Residuals 


o  Expected  Value  (average  reeidual)  ...........  0.294  FT 

o  Standard  Deviation  . .  5 . 934  FT 

o  Average  Abaolute  Error  . . .  3 . 595  FT 

o  Maximum  Abaolute  Error . .  99.610  FT 

o  Percent  Absolute  Values  <  12.5  FT  . . 95.141  Z 

o  Number  Absolute  Values  >99.5  FT  ...........  1 

o  Number  of 


Positive  Residuals  (overestimates)  ...  222  417 
Negative  Reslduals(underestimates)  . . .  205  365 


Zero  Residuals  . . .  38  678 

*> 

o  Maximum  Absolute  Adjusted  Residual  .  65.588  FT  - 

o  Average  Absolute  Adjusted  Residual . .  3.12-  FT 


o  Percent  Absolute  Adjusted  Values  >  40  FT  ....  0.020  Z 

The  adjusted  residuals  are  seen  to  be  much  smaller  than  the  unadjusted  ones. 
This  indicates  that  big  values  of  the  (unadjusted)  residuals  are  largely 
confined  to  steep  terrain.  Results  for  the  second  run  are  displayed  below.  It 
is  seen  that  not  to  use  tangential  information  results  in  a  definite 
deterioration  of  accuracy: 


o  Expected  Value  (average  residual)  .  0.209  FT 

o  Standard  Deviation  . .  6.034  FT 

o  Average  Absolute  Error  . 3.657  FT 

o  Maximum  Absolute  Error  . .  160.619  FT 

o  Percent  Absolute  Values  <  12.5  FT  . .  95.044  Z 

o  Number  Absolute  Values  >99.5  FT  .  15 

o  Number  of 


Positive  Residuals  (overestimates)  ...  221  420 
Negative  Residuals(underestimates)  ...  206  419 
Zero  Residuals  . . .  38  621 
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As  expected,  the  third  run,  featuring  linear  interpolation,  did  not  achieve  the 
accuracy  of  the  first  run.  It  comes  as  somewhat  cf  a  surprise,  however,  that 
the  loss  of  accuracy  is  not  more  pronounced.  An  analysis  of  the  results  shows 
that  this  is  due  to  the  extraordinary  large  number  of  zero  residuals.  This 
phenomenon,  in  turn,  is  an  artefact  of  the  digitization:  digitized  contour  lines 
contain  many  groups  of  successive  points  that  lie  on  a  common  line.  If  two  of 
such  points  are  vertices  of  the  same  triangle,  then  they  and  all  intermediate 
points  reproduce  elevation  under  linear  interpolation.  As  we  pointed  out 
earlier,  the  precise  reproduction  of  digitized  points  is  not  necessarily 
desirable  because  the  latter  carry  digitization  errror.  Below  are  the 
statistics  for  the  linear  run: 


o  Expected  Value  (average  residual)  .  0.326  FT 

c  Standard  Deviation  . . .  5.856  FT 

o  Average  Absolute  Error  . .  3.557  FT 

o  Maximum  Absolute  Error  . 106.557  FT 

o  Percent  Absolute  Values  <  12.5  FT  .  95.277  % 

o  Number  Absolute  Values  >  99.5  FT  .  4 

o  Number  of 


Positive  Residuals  (overestimates)  ...  190  778 
Negative  Residuals(underestimates)  ...  175  542 
Zero  Residuals  . .  100  140 

6 . 3  CONCLUSIONS 


Two  particular  concerns  were  our  capability  of  obtaining  the  Voronoi 
triangulation  of  a  sufficiently  large  set  of  sample  points  and  of  solving  the 
large  linear  system  for  the  gradients  at  the  sample  points,  using  commonly 
available  computer  resources  With  respect  to  these  concerns  we  found  that 
samples  of  up  to  70,000  points  could  be  triangulated  within  a  reasonable  time 
frame.  In  fact,  CPU  time  has  been  less  of  a  problem  than  memory  space,  a 
limitation  that  was  overcome  by  developing  a  decomposition  method.  The 
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Gauss-Seidel  method  described  in  Section  5.2  was  found  to  be  faster  than 
anticipated,  both  with  respect  to  the  CPU  time  required  by  a  single  iteration 
(approximately  150  CPU  seconds/iteration)  and  the  speed  of  convergence  (ten 
Iterations).  Alternative  methods  will  perform  as  well  or  faster. 

To  sum  up,  we  conclude  that  the  computational  effort  of  using  nonlinear 
techniques  for  the  generation  of  a  smooth  synthetic  surface  and  subsequent 
regular  grid  is  substantial,  but  within  the  bounds  for  routine  calculations  on  a 
computer  of  medium  size  such  as  the  VAX  11/780.  Thus  we  were  succesful  in 
achieving  the  major  goals  of  this  feasibility  study. 

The  residuals  obtained  when  trying  to  recover  the  original  digitized  contour 
points  are  by  and  large  comparable  to  the  resolution  of  the  digitized  data. 
However,  there  are  instances  of  very  large  residuals  which  may  well  be 
unacceptable  for  subsequent  applications.  Such  discrepancies  were  expected 
because  our  method  at  this  point  still  lacks  the  capability  to  handle 
cartographic  features  at  which  the  actual  terrain  surface  is  not  smooth,  that 
is,  it  exhibits  discontinuities  of  slope.  Such  features  include  lake  shores, 
river  banks,  as  well  as  some  ridge  and  drainage  lines.  Modeling  by  a  smooth 
surface  may  not  be  sufficiently  accurate  at  such  locations. 

In  order  to  achieve  the  full  accuracy  of  our  approach,  the  following  measures 
need  to  be  taken: 

o  Cancel  the  smoothness  requirements  along  lake  shores,  river  banks,  as 
well  as  along  certain  ridge  and  drainage  lines. 

o  Smooth  digitized  contour  lines  and  determine  tangent  directions  prior 
to  thinning. 

o  Investigate  adaptive  tolerance  selection  schemes  for  higher  density 
sampling  in  rough  terrain. 


o  Compare  local  and  global  methods  for  specifying  tangent  planes 
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